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1.0  INTRODUCTION 

Rational  design  and  advancement  in  materials  science  will  ultimately  rely  on  an  atomic-scale 
understanding  of  the  targeted  functionality.  Corresponding  modeling  must  then  address  the  behavior  of 
electrons  and  the  resulting  interactions  (often  expressed  in  the  terminology  of  chemical  bonds)  that  govern 
the  elementary  processes  among  the  atoms  and  molecules  in  the  system.  Modem  electronic  structure 
theory  methods  like  density-functional  theory  (DFT)  [1-5]  have  matured  to  a  standard  tool  for  this  task, 
allowing  a  description  that  is  often  already  accurate  enough  to  allow  for  a  modeling  with  predictive 
character.  These  techniques  are  referred  to  as  first-principles  (or  in  latin:  ab  initio )  to  indicate  that  they  do 
not  rely  on  empirical  or  fitted  parameters,  which  then  makes  them  applicable  for  a  wide  range  of  realistic 
conditions,  e.g.  realistic  environmental  situations  of  varying  temperatures  and  pressures  [6]. 

The  latter  type  of  application  seems  at  first  sight  at  variance  with  the  frequent  argument  describing  DFT  as 
a  zero-temperature,  zero-pressure  technique.  Such  a  confusion  arises,  when  thinking  that  DFT  provides 
(apart  from  a  wealth  of  information  about  the  electronic  structure)  “only”  the  total  energy  of  the  system. 
Instead,  it  is  crucial  to  realize  that  this  kind  of  energetic  information  can  be  obtained  as  a  function  of  the 
atomic  configuration  {Ry}.  This  leads  to  the  so-called  potential  energy  surface  (PES)  ii({Ry}),  which  then 
contains  the  relevant  information  needed  to  describe  the  effect  of  temperature  on  the  atomic  positions. 
Obviously,  a  (meta)stable  atomic  configuration  corresponds  to  a  (local)  minimum  of  the  PES.  The  forces 
acting  on  the  given  atomic  configuration  are  just  the  local  gradient  of  the  PES,  and  the  vibrational  modes 
of  a  (local)  minimum  are  given  by  the  local  PES  curvature  around  it. 

One  possibility  to  go  from  this  to  situations  of  finite  temperature  and  finite  pressure  is  to  achieve  a 
matching  with  thermodynamics.  This  is  the  general  idea  behind  ab  initio  atomistic  thermodynamics, 
namely  to  employ  the  information  on  the  first-principles  PES  to  calculate  appropriate  thermodynamic 
potential  functions  like  the  Gibbs  free  energy  [7-10].  Once  such  a  quantity  is  known,  one  is  immediately  in 
a  position  to  evaluate  macroscopic  system  properties  using  the  standard  methodology  of  thermodynamics. 
Apart  from  bridging  to  any  (T/fi-conditions,  this  methodology  is  particularly  useful  for  larger  systems, 
which  may  readily  be  divided  into  smaller  subsystems  that  are  mutually  in  equilibrium  with  each  other. 
Each  of  the  smaller  and  thus  potentially  simpler  subsystems  can  then  be  first  treated  separately,  and  the 
contact  between  the  subsystems  is  thereafter  established  by  relating  their  corresponding  thermodynamic 
potentials.  Such  a  “divide  and  conquer”  type  of  approach  can  be  especially  efficient,  if  infinite,  but 
homogeneous  parts  of  the  system  like  bulk  or  surrounding  gas  phase  can  be  separated  off,  and  are  then 
merely  represented  by  corresponding  reservoirs  [11-16].  Although  not  further  discussed  here,  another 
aspect  could  be  to  consider  situations  of  “constrained  equilibria”  [14,15],  where  not  all,  but  only  some  of 
the  subsystems  are  in  thermodynamic  equilibrium. 
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Figure  1:  Schematic  representation  of  the  system  discussed  here,  a  single-crystal  metal  surface 
in  contact  with  a  surrounding  gas  phase  characterized  by  defined  temperature  T  and  pressure  p. 
The  shaded  area  represents  the  finite  part  of  the  system  that  is  affected  by  the  presence  of  the 
surface. 


2.0  SURFACES  IN  REALISTIC  ENVIRONMENTS 

2.1  Surface  free  energy 

Here  we  will  illustrate  how  this  quite  general  concept  works  and  what  it  can  contribute  in  practice  by 
using  it  to  determine  the  equilibrium  geometry  and  composition  of  a  solid  surface  in  contact  with  a  given 
environment  at  finite  temperature  and  pressure.  For  the  sake  of  simplicity  we  consider  the  case  of  a 
monoatomic  metal  and  an  oxygen  atmosphere  described  by  an  oxygen  pressure  p  and  a  temperature  T  [17- 
20],  and  refer  to  the  literature  for  the  generalizations  to  compounds  like  oxides  [13]  or  alloys  [21],  and  to 
environments  containing  multiple  gas  phase  species  [14,15].  Under  conditions  of  defined  ( T,p ),  the 
appropriate  thermodynamic  potential  to  consider  is  the  Gibbs  free  energy.  For  this  quantity  we  introduce 
the  following  notation:  Capital  G  refers  to  an  absolute  Gibbs  free  energy,  while  lower  case  g  is  used  to 
denote  a  Gibbs  free  energy  per  formula  unit  or  particle.  In  the  case  of  an  infinite,  homogeneous  system  the 
latter  is  then  equivalent  to  the  chemical  potential  p,  i.e.  if  the  homogeneous  system  is  viewed  as  a 
reservoir,  p  gives  the  cost  at  which  this  reservoir  provides  particles.  For  convenience  we  will  use  the 
symbol  p  instead  of  g,  when  we  want  to  emphasize  a  system’s  role  as  a  reservoir.  In  the  present  context  it 
is  particularly  the  oxygen  environment  which  acts  as  such  a  reservoir,  because  it  can  give  (or  take)  any 
amount  of  oxygen  to  (or  from)  the  sample  without  changing  the  temperature  or  pressure. 

Figure  1  shows  a  schematic  representation  of  the  system  discussed  here,  a  solid  phase  in  contact  with  a 
surrounding  gas  phase.  We  can  break  down  the  Gibbs  free  energy  of  this  entire  system  into  contributions 
coming  from  the  bulk  of  the  solid  phase  Gsoiid,  from  the  homogeneous  gas  phase  Ggas  and  an  additional 
term  introduced  through  the  surface  AGsurf 

(1)  G  Gsolid  Ggas  AGsurf 
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If  the  surface  is  homogeneous  as  in  the  case  of  an  ideal  single-crystal  surface,  AGsurf  will  scale  linearly 
with  the  surface  area  A,  and  we  can  introduce  the  surface  free  energy  per  unit  area,  y.  Upon  rearranging  eq. 

(1)  we  therefore  have 

(2)  y  =  —  ( G  -  Gsolid  -  Gms  ) 

A 

We  notice  that  y  is  well  defined  by  a  finite  part  of  the  total  (infinite)  system.  With  increasing  distance 
from  the  surface,  eventually  both  the  solid  and  the  gas  phase  part  of  the  total  system  will  no  longer  be 
affected  by  the  created  surface.  Although  contained  in  G,  these  (infinite)  parts  of  the  total  system  are  then 
effectively  canceled  out  in  eq.  (2)  by  the  subtraction  of  the  equivalent  amounts  of  homogeneous  systems 
(GgaS  and  GSOHd).  We  can  therefore  concentrate  on  the  finite  part  of  the  system  that  is  affected  by  the 
surface.  If  this  part  contains  Am  metal  atoms  and  No  oxygen  atoms  per  surface  area,  this  allows  us  to 
rewrite  eq.  (2)  as 

(3)  y{T,p)  =  UG(T,p,N0,Nu)-Nugu{T,p)-N0fi0{T,p ))  , 

A 

where  we  have  introduced  the  Gibbs  free  energy  per  metal  atom  gM  in  the  bulk,  and  the  oxygen  chemical 
potential  u()  of  the  gas  phase. 

At  this  stage  it  is  appropriate  to  spend  a  few  words  on  the  sign  convention.  In  this  text,  a  more  negative 
Gibbs  free  energy  will  indicate  a  more  stable  state  of  the  system.  In  the  interpretation  of  a  chemical 
potential  this  translates  to  p  approaching  -go  in  the  limit  of  an  infinitely  dilute  gas,  since  adding  a  particle 
will  then  yield  an  infinite  gain  in  entropy.  As  a  consequence,  y  >  0  indicates  the  cost  of  creating  the 
surface  between  the  solid  bulk  phase  and  the  homogeneous  gas  phase.  Alternatively,  when  discussing  the 
stability  of  phases  that  result  from  adsorbing  species  at  the  solid  surface,  it  can  be  convenient  to  choose 
another  zero  reference.  Instead  of  y  as  a  measure  of  the  cost  of  creating  the  surface  at  all,  one  can 
introduce  the  surface  free  energy  of  the  clean  surface 

(4)  yaew(T,p)  =  4  (  G(T ,p,0,NM)-  NMgM (T,p)) 

A 

and  evaluate  the  Gibbs  free  energy  of  adsoiption  as  a  measure  of  the  cost  with  respect  to  the  clean  surface 
[17,19,20] 

A  G*d(T,p)  =  yclew(T,P,0,N'M)-y(T  ,p,N0,NM)  = 

(5)  =  (G(T,p,No,NM)-G(T,p,0,N'M)-Nopo(T,p)-{NM-N'M)gM{T,p)) 

where  the  last  term  accounts  for  a  possible  difference  in  the  number  of  metals  atoms  between  the  reference 
clean  surface  and  the  oxidized  surface  structural  model.  Obviously,  the  most  stable  surface  structure  and 
composition  at  given  (T,p)  in  the  gas  phase  is  the  one  that  minimizes  the  surface  free  energy,  or 
equivalently  the  one  that  leads  to  the  most  positive  Gibbs  free  energy  of  adsorption  at  the  corresponding 
oxygen  chemical  potential. 

2.2  Calculating  Gibbs  free  energies 

As  apparent  from  eqs.  (3)  and  (5)  the  quantities  determining  y  or  AGad  are  the  Gibbs  free  energy  of  the 
solid  surface  and  of  the  solid  bulk,  as  well  as  the  chemical  potential  of  the  oxygen  environment.  Since  the 
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contributions  to  free  energies  in  the  gas  phase  and  in  the  solid  phase  are  quite  different,  we  will  discuss  the 
evaluation  of  these  quantities  separately.  First  we  concentrate  on  the  oxygen  chemical  potential,  and  then 
on  the  computation  of  the  solid  phase  Gibbs  free  energies. 


2.2.1  Gas  phase  chemical  potential 

The  chemical  potential  of  oxygen,  which  enters  eqs.  (3)  and  (5),  is  determined  by  the  condition  of 
thermodynamic  equilibrium  with  the  surrounding  gas  phase  reservoir.  At  the  accuracy  level  relevant  for 
our  purposes,  this  gas  is  well  described  by  ideal  gas  laws.  If  we  write  the  chemical  potential  at  given 
temperature  T  and  pressure  p  as 

(6)  fiQ (T, p)  =  y2  Ho2 (gas) (T,p)  =  y2(-kBT In Qqi(bu)  +  pV)i N 

we  therefore  have  to  evaluate  the  partition  function  0o*(gaS)  of  an  ideal  gas  composed  of  N 
indistinguishable  CF  molecules 

(7)  0o2(gas)  “l?oJ  —l?  q  q  q  q  ) 


with  q0  the  partition  function  of  one  CF  molecule,  which  can  further  be  subdivided  into  different 

partition  functions  with  obvious  nomenclature.  In  writing  eq.  (7)  we  assumed  that  the  nuclear/electronic 
(nucl,  electr)  degrees  of  freedom  are  decoupled  from  the  v  i  b  r  a  t  i  o  n  a  I  /  r  o  t  a  t  i  o  n  a  1  (vib,  rot)  ones  (Bom- 
Oppenheimer  approximation),  and  further  that  also  vibrational  and  rotational  motions  are  decoupled  as 
they  take  place  on  different  time  scales. 

Inserting  eq.  (7)  into  eq.  (6)  we  then  arrive  at 


(8)  p0(T,p)  =  -— 


(  i  Y 

—(qtmm)N 

kBT  In 

-pV 

{N!  ) 

+  Yi  vrot  +  X  A  +  X  A  +'Am 


nucl 


and  proceed  here  by  briefly  commenting  on  the  translational,  rotational,  vibrational,  electronic  and  nuclear 
free  energy  terms.  A  more  detailed  derivation  can  be  found  in  most  textbooks  on  Statistical  Mechanics, 
e.g.  the  one  by  Me  Quarrie  [22]. 


Translational  free  energy 

In  the  classical  limit  the  energy  due  to  the  center-of-mass  motion  of  a  particle  in  a  box  is 

?;2k2  .  n  (  „  „  „\ 

(9)  £k  = -  ,  k  =  —  In Xx  +  nvy  +  nzz] , 

2m  Ly  •  ’ 

where  L  =  Vin  characterizes  the  box  size,  m  is  the  mass  of  one  particle,  x/  y  /  z  are  unit- vectors  in  the 

three  cartesian  directions,  and  «x/y/z  go  from  1  to  go.  In  the  thermodynamic  limit  ( L  — mx>),  the  one-particle 
partition  function  becomes 

f  2nmkBT  V/2 


(10)  trans  _  y 
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With  this  we  rewrite  the  term  in  brackets  in  eq.  (8)  to  (employing  the  ideal  gas  law  pV  =  NkBT  and  the 
Stirling  formula  In  N\  ~  N  In/V  -  N  at  some  stage) 


(11) - 

2N 


kBT  In 

(  1  ) 

1  (^tranS)^ 

-pV 

=  -y2kBT  In 

z  2n m N 

3/2  ( kBT f2 

- 

l  N\  ) 

V  h2  ) 

P 

Rotational  free  energy 

In  the  rigid  rotator  approximation  the  rotational  partition  function  is  written  as 


(12)  qm  = 


j= o 


where  B0  =  h2/(2I)  is  the  rotational  constant  and  /  the  moment  of  inertia  depending  primarily  on  mass  and 
equilibrium  bond  lengths.  In  case  of  homonuclear  diatomic  molecules  like  CF  (or  other  molecules  that 
have  multiple  indistinguishable  orientations),  qmt  is  a  bit  more  tricky  as  it  couples  with  the  nuclear  spin 
degrees  of  freedom  (the  total  wave  function  must  be  antisymmetric  under  exchange  of  the 
indistinguishable  particles).  At  the  temperatures  of  interest  to  us,  this  can  be  approximately  combined  into 
a  classical  symmetry  number  <rsym  indicating  the  number  of  indistinguishable  orientations  that  the  molecule 
can  have  (e.g.  =1  for  heteronuclear  diatomic  molecules,  =2  for  homonuclear  diatomic  molecules).  At  such 
temperatures,  where  the  spacing  of  the  rotational  levels  is  small  compared  to  kBT,  the  sum  in  qmt  can  be 
converted  into  an  integral  with  the  Euler-Maclaurin  series  and  one  ends  up  with 


(13) 


plot  «  ~kBT  In 


ytjsymB0  j 


Notice  that  this  holds  only  for  linear  molecules,  where  then  only  B0  enters.  In  more  complex  cases  one 
would  need  to  diagonalize  the  inertial  tensor  and  consider  all  three  eigenvalues  Aa,  B0  and  C0. 


Vibrational  free  energy > 

The  vibrational  contribution  is  obtained  within  the  harmonic  approximation  by  writing  the  partition 
function  as  a  sum  over  the  harmonic  oscillators  of  all  M  fundamental  modes  a>\  of  the  particle, 


(14) 


„  vib 


=  ZZexp 


(=1  72=0 


-  (»  + 

kBT 


Evaluation  of  the  geometric  series  yields 

M  f- 

(15)  //ib=£ZPE+A//lb  =  £ 

!=i 


kB7Tn 


1  -  exp 


f t  W 
nco 


V^b  T  JJ 


where  the  first  term  arises  from  the  zero-point  vibrations. 

Electronic  and  nuclear  free  energy 

For  most  molecules  internal  excitation  energies  are  large  compared  to  kBT,  so  that  the  only  term 
contributing  significantly  to  the  partition  function  is  the  ground  state.  Taking  a  possible  spin  degeneracy  of 
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this  ground  state  into  account,  we  end  up  with 
( 1 6)  //lectr  «  £“al  -  kBT In (/spin ) 

i.e.  with  the  total  energy  of  the  02  molecule  and  a  term  depending  on  the  electronic  spin  degeneracy  of  the 
ground  state,  /pm.  The  same  form  would  also  be  obtained  for  the  nuclear  degrees  of  freedom,  which  are 
even  more  confined  to  the  ground  state  due  to  the  much  larger  separation  of  nuclear  energy  levels. 
Although  the  nuclear  partition  function  may  thus  be  different  from  unity  (according  to  the  degeneracy  of 
the  nuclear  spin  ground  state),  we  will  omit  it  here,  since  the  nuclear  state  is  rarely  altered  in  chemical 
processes  and  therefore  does  not  contribute  to  the  thermodynamic  changes  discussed  here. 


Figure  2:  Temperature  dependence  of  the  relative  oxygen  chemical  potential  Ajuo(T,1  atm)  at  p=1 
atm.  Compared  are  tabulated  values  from  ref.  [23]  (crosses)  with  the  result  of  eq.  (18)  using  the 
material  parameters:  orSym=2,  /spm=3,  Bo=0.18  meV,  and  cu0=196  meV  [23].  Additionally  shown  are 
the  sums  of  the  individual  contributions:  vibrational  (dashed  line,  almost  coinciding  with  the 
zero  axis),  vibrational+nuclear  (dotted  line),  vibrational+nuclear+rotational  (dash-dotted  line). 
The  remaining  large  difference  to  the  full  result  (solid  line)  is  due  to  the  translational 
contribution. 


Bringing  it  all  together 

Equations  (11),  (13),  (15),  (16)  allow  the  analytic  evaluation  of  the  terms  entering  eq.  (8),  based  on  the 
total  energy  and  vibrational  modes  of  the  gas  phase  molecule,  as  well  as  the  molecular  mass,  the  rotational 
constant,  the  symmetry  number  and  the  spin  degeneracy  (which  are  available  in  thermochemical  tables 
[23]).  As  will  become  clear  below,  it  is  convenient  to  separate  out  the  total  and  zero  point  energy  terms 
and  write  eq.  (8)  in  the  following  form 

(17)  n0  (T ,  p)  =  y2E^l  +  y2E^  +  A Mo(T,p) 

where  A juo(T,p)  contains  now  all  temperature  and  pressure  dependent  free  energy  contributions 
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(18)  Ap0  (T,p)  =  -y2kBT\\n 


2mn 


\  3/2 


(kBT) 


5/2 


+  In 


f  kBT  ^ 
V*symB0; 


-In 


1  -  exp 


^  Tico.  ^ 


■  ln(/spin)l 


From  the  structure  of  eq.  ( 1 8)  it  is  apparent  that  one  can  rewrite  eq.  ( 1 7)  in  the  following  form 


(19) 


Mo(T,p)  =  y2ESm+y2E™E+AMc 


(T,p°)  +  ykBT\n 


„o 

\  IJ  y 


which  allows  to  use  tabulated  enthalpy  and  entropy  values  at  standard  pressure  p°=  latm  [23]  to  determine 
A po(T,p)  instead  of  using  eq.  (18)  [13].  For  oxygen,  both  approaches  yield  virtually  identical  results  in  the 
temperature  range  discussed  here,  as  illustrated  in  Fig.  2  for  Ap0(T,  1  atm).  However,  for  other  (more 
complex)  molecules  one  or  the  other  approach  may  be  more  convenient,  depending  on  the  availability  of 
the  gas  phase  data. 


2.2.2  Gibbs  free  energy  of  solid  bulk  and  surface 

Similar  to  the  procedure  applied  for  the  gas  phase  chemical  potential  we  address  the  computation  of  the 
solid  phase  Gibbs  free  energies  by  first  decomposing  them  into  several  contributing  terms  [13,15] 

(20)  G  =  £total  +  F’vib  +  Econf+ p  V 

where  Z/llllal  is  the  total  (internal)  energy,  F"b  the  vibrational  free  energy,  and  A0"1  the  configurational  free 
energy.  A  crucial  aspect  that  governs  our  analysis  of  all  of  these  terms  is  that  the  quantities  of  interest  to 
us,  namely  surface  free  energies  or  Gibbs  free  energies  of  adsorption,  do  not  depend  on  absolute  Gibbs 
free  energies.  What  enters  into  the  corresponding  eqs.  (3)  or  (5)  are  differences  of  Gibbs  free  energies 
only.  This  can  allow  for  quite  some  degree  of  error  cancellation,  in  particular  if  the  different  terms 
correspond  to  rather  similar  situations  like  the  Gibbs  free  energies  of  solid  bulk  and  solid  surface. 

The  dominant  term  in  eq.  (20)  is  the  total  energy,  which  as  discussed  in  the  introduction  is  provided 
through  the  DFT  calculations.  We  note  in  passing  that  the  thermodynamic  formalism  is,  of  course,  general 
and  would  be  equally  valid  when  using  total  energies  coming  from  e.g.  less  accurate  schemes.  However, 
when  using  say  a  semi-empirical  potential  to  provide  the  total  energies,  the  accuracy  of  the  ensuing 
thermodynamic  reasoning  would  also  only  be  at  this  level.  It  is  precisely  the  idea  of  ab  initio  atomistic 
thermodynamics  to  carefully  evaluate  the  total  energy  contributions  and  thereby  carry  over  the  predictive 
power  of  the  first-principles  technique  to  finite  (Tp)-conditions. 

Turning  to  the  other  terms  in  eq.  (20),  we  find  from  a  simple  dimensional  analysis  that  its  contribution  to 
the  surface  free  energy  (normalized  to  unit  surface  area)  will  be  [pV/A]  =  atm  A3/A2  ~  1 0  3  meV/A2.  Even 
for  p  ~  100  atm,  the // E-contribution  will  therefore  still  be  less  than  ~0.1  meV/  A2.  We  will  see  below  that 
in  the  intended  application  to  metal  surfaces  in  contact  with  realistic  environments  this  is  a  negligible 
contribution  compared  to  the  other  free  energy  terms,  which  are  of  the  order  of  tens  of  meV/  A2. 

The  contribution  from  the  vibrational  degrees  of  freedom  can  be  handled  with  the  same  harmonic 
approximation  applied  already  for  the  gas  phase  chemical  potential.  However,  instead  of  writing  the 
vibrational  free  energy  as  arising  from  a  sum  of  discrete  fundamental  modes  to,-,  cf.  eq.  (15),  we  now 
introduce  the  phonon  density  of  states  (DOS)  c(co)  and  obtain 

(21)  Fvib  =  | dco  Fvib  (T,a>)  a(co) 
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Figure  3:  Estimated  vibrational  contribution  AFvlb'ad  to  the  Gibbs  free  energy  of  adsorption  for 
the  p( 2x2)  O/Pd(100)  structure.  Used  is  the  Einstein  approximation  to  the  phonon  density  of 
states,  where  the  characteristic  frequency  of  O  atoms  is  changed  from  196  meV  in  the  gas  phase 
to  40  meV  (black  line),  80  meV  (blue  line)  or  120  meV  (red  line)  at  the  surface.  Changes  of  the 
vibrational  modes  of  surface  Pd  atoms  upon  O  adsorption  are  neglected. 
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A  proper  evaluation  of  the  vibrational  contribution  to  AGad  or  y  amounts  therefore  to  computing  the 
phonon  DOS  of  the  solid  bulk  and  surface.  This  information  is  contained  in  the  PES  and  can 
correspondingly  be  calculated  by  DFT  [5].  Elowever,  since  surface  phonon  DOS  calculations  are 
computationally  still  quite  involved,  getting  away  with  simpler  approximations  would  be  particularly 
worthwhile.  For  this,  we  come  back  to  the  afore  mentioned  cancellation  in  the  differences  of  Gibbs  free 
energies.  It  is  crucial  to  realize  that  the  vibrational  free  energy  of  the  solid  bulk  and  surface  is  not  small  or 
negligible,  yet  what  matters  for  the  determination  of  e.g.  the  Gibbs  free  energy  of  adsorption  is  only  the 
following  difference 


Apvib.ad^ 
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where  a  is  the  phonon  DOS  of  the  surface  considered,  <7ciCan  the  one  of  the  corresponding  clean  surface,  <7m 
the  one  of  the  metal  bulk,  and  in  the  second  line  we  have  realized  that  the  temperature  dependent 

vibrational  contribution  to  piQ^(T)  is  negligible,  cf.  Fig.  2.  Equation  (23)  shows  that  only  the  changes  in 

the  vibrational  properties  of  the  atoms  at  the  surface  compared  to  their  counterparts  in  the  clean  surface, 
the  solid  bulk  or  the  gas  phase  are  decisive.  Before  initiating  a  full-blown  surface  phonon  DOS 


2-8 


RTO-EN-AVT-142 


Ab  Initio  Atomistic  Thermodynamics  for  Surfaces:  A  Primer 


calculation,  it  can  therefore  be  very  valuable  to  use  approximations  to  a,  adcd„,  and  crM ,  and  obtain  an  order 
of  magnitude  estimate  of  AFv,b'dd  first  [13].  A  very  simple  approximation  would  e.g.  be  an  Einstein  model, 
which  considers  only  one  characteristic  frequency  for  each  atom  type.  Allowing  this  frequency  to  change 
significantly  for  surface  atoms  compared  to  those  in  the  solid  bulk  or  in  the  gas  phase  provides  then  a  first 
idea  of  the  magnitude  of  the  vibrational  contribution.  This  approach  is  exemplified  in  Fig.  3  for  the  below 
discussed  case  of  a  Pd(100)  surface  in  contact  with  an  oxygen  enviromnent.  Even  for  quite  extensive 
variations  of  the  characteristic  vibrational  modes,  the  resulting  AF"  ,ad  always  stays  within  ~  ±5  meV/  A2 
for  the  entire  temperature  range  up  to  600  K  that  is  of  interest  to  our  study.  Comparing  with  the  much 
larger  contribution  from  the  total  energy  terms,  we  will  use  this  information  to  justify  neglecting  the 
vibrational  AFvlbad  below.  However,  we  stress  that  this  is  not  a  general  result.  There  might  well  be 
applications  where  the  inclusion  of  vibrational  effects  on  the  surface  free  energy  or  Gibbs  free  energy  of 
adsorption  can  be  important  as  e.g.  in  the  adsorption  of  larger  molecules  [24]  or  in  hydrogen  containing 
environments  [25].  In  such  cases,  it  may  already  be  sufficient  to  only  consider  some  prominent  vibrational 
modes,  but  eventually  an  explicit  calculation  of  the  surface  and  bulk  phonon  DOS  may  be  required. 

This  leaves  as  remaining  term  the  configurational  free  energy.  A  full  evaluation  of  this  contribution  is 
computationally  most  involved,  since  it  requires  a  proper  sampling  of  the  huge  configuration  space 
spanned  by  all  possible  surface  structures.  Although  modem  statistical  mechanics  methods  like  Monte 
Carlo  simulations  [26,27]  are  particularly  designed  to  efficiently  fulfill  this  purpose,  they  still  require  a 
prohibitively  large  number  of  free  energy  evaluations  to  be  directly  linked  with  electronic  structure 
theories  [6],  A  way  to  circumvent  this  problem  is  to  map  the  real  system  somehow  onto  a  simpler, 
typically  discretized  model  system,  the  Hamiltonian  of  which  is  sufficiently  fast  to  evaluate.  Obvious 
uncertainties  of  this  approach  are  how  appropriate  the  model  system  represents  the  real  system,  and  how 
its  parameters  can  be  determined  from  the  first-principles  calculations.  The  advantage,  on  the  other  hand, 
is  that  such  a  detour  via  an  appropriate  (“coarse-grained”)  model  system  often  provides  deeper  insight  and 
understanding  of  the  ruling  mechanisms.  If  the  considered  problem  can  be  described  by  a  lattice  defining 
e.g.  the  possible  adsorption  sites  for  the  gas  phase  species  in  the  system,  a  prominent  example  for  such  a 
mapping  approach  comes  under  the  names  lattice-gas  Hamiltonians  (LGHs)  or  cluster  expansions  (CEs) 
[6,28-30]. 

Here,  we  will  instead  concentrate  on  a  much  simpler  alternative,  which  focuses  on  screening  a  number  of 
known  (or  possibly  relevant)  ordered  surface  structures  by  directly  comparing  which  of  them  turns  out  to 
be  most  stable  under  which  (7’, /^-conditions,  i.e.  which  of  them  exhibits  the  lowest  surface  free  energy  or 
Gibbs  free  energy  of  adsorption.  For  sufficiently  low  temperatures,  the  remaining  configurational  entropy 
per  surface  area  is  then  only  due  to  a  limited  number  of  defects  in  these  ordered  structures  and  can  be 
estimated  as 


(24) 
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where  n  is  the  small  number  of  defects  in  a  system  with  N  surface  sites  (n  «  TV).  For  N,  n  »  1  we  can 
apply  the  Stirling  formula  which  leads  to 
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With  a  typical  area  per  surface  site  of  A  ~  10  A2  for  transition  metal  surfaces,  we  correspondingly  deduce 
that  the  configurational  entropy  contribution  to  the  Gibbs  free  energy  is  less  than  3  meV/  A2  for  any  T  < 
1000  K  [15].  We  will  see  below  that  in  the  application  to  a  Pd  surface  in  contact  with  an  oxygen 
environment  this  is  almost  always  negligible  and  will  qualify  to  which  changes  its  explicit  consideration 
will  lead.  While  the  effect  of  configurational  entropy  due  to  disorder  is  therefore  not  of  much  concern,  the 
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obvious  limitation  in  this  direct  screening  approach  is  that  its  reliability  is  restricted  to  the  number  of 
considered  configurations,  or  in  other  words  that  only  the  stability  of  those  structures  plugged  in  can  be 
compared.  The  predictive  power  extends  therefore  only  to  those  structures  that  are  directly  considered,  i.e. 
the  existence  of  unanticipated  surface  geometries  or  stoichiometries  cannot  be  predicted.  As  such, 
appropriate  care  should  be  in  place  when  addressing  systems  where  only  limited  information  about  the 
surface  structures  is  available.  With  this  in  mind,  the  direct  screening  approach  to  ab  initio  atomistic 
thermodynamics  can  still  be  a  particularly  valuable  tool,  since  it  allows,  for  example,  to  rapidly  compare 
the  stability  of  newly  devised  structural  models  against  existing  ones.  In  this  way,  it  gives  tutorial  insight 
into  what  structural  motives  may  be  particularly  important,  which  may  even  yield  ideas  about  other 
structures  that  one  should  test  as  well.  Still,  the  limited  reliability  to  the  set  of  actually  considered 
structural  models  must  always  be  borne  in  mind  and  can  really  only  be  overcome  by  a  proper  sampling  of 
configurational  space,  which  then  leads  also  to  a  more  general  and  systematic  way  of  treating  phase 
coexistence  and  order-disorder  transitions. 
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Figure  4:  Generic  free  energy  plot  (a)  and  surface  phase  diagram  (b)  for  a  surface  in  equilibrium 
with  a  surrounding  oxygen  gas  phase,  a)  An  adsorbate  phase  will  become  more  stable  than  the 
zero  reference  clean  metal  surface,  if  its  AGad  >0  for  some  oxygen  chemical  potential  (note  the 
inverted  y-scale!).  If  there  is  more  than  one  adsorbate  phase,  always  the  one  with  the  largest 
AGad  will  become  most  stable,  as  indicated  here  by  the  red  line.  Finally,  for  Apo  >  1/y  A//^r=0K), 
the  bulk  oxide  will  always  result  as  the  most  stable  phase,  b)  Converting  the  obtained  A p0 
stability  range  for  each  phase  (indicated  schematically  by  the  different  shaded  regions  in  the 
free  energy  plot  on  the  left)  into  (T,p)-conditions  using  eq.  (18)  allows  to  draw  the  resulting 
surface  phase  diagram. 
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2.3  Free  energy  plots  and  surface  phase  diagrams 

We  will  now  employ  the  general  thermodynamic  framework  developed  in  sections  2.1  and  2.2  to 
investigate  the  structure  and  composition  of  a  solid  surface  in  contact  with  a  given  environment  at  finite 
temperature  and  pressure.  More  specifically,  we  aim  at  applying  it  to  describe  a  single-crystal  transition 
metal  surface  in  contact  with  an  oxygen  gas  phase,  but  before  we  address  the  actual  case  of  a  Pd(100) 
surface  below  let  us  first  spend  some  time  with  more  general  considerations.  In  light  of  the  discussion  in 
section  2.2.2  we  will  neglect  the  vibrational  contribution  A/-’v,b'ad  to  the  Gibbs  free  energy  of  adsorption, 
and  within  the  spirit  of  the  direct  screening  approach  we  will  also  neglect  the  configurational  entropy  term 
in  the  solid  phase  Gibbs  free  energies  for  the  time  being.  This  transforms  the  general  eq.  (5)  into  the 
working  equation 


(26)  \Gad(T,p) « 
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Evaluating  A Gad  requires  now  primarily  total  energies  which  are  directly  amenable  to  electronic  structure 
theory  calculations.  In  the  direct  screening  approach  to  ab  initio  atomistic  thermodynamics,  these 
quantities  would  be  calculated  for  a  number  of  ordered  surface  structural  models.  Equation  (26)  allows 
then  to  directly  plot  the  Gibbs  free  energy  of  adsorption  for  each  model  as  a  function  of  the  oxygen 
chemical  potential  A p0,  as  illustrated  schematically  in  Fig.  4a.  This  yields  a  straight  line  for  each  model 
considered,  and  at  any  given  A p0  the  model  with  the  lowest  lying  line  (most  stable  AGad)  is  identified  as 
the  most  stable  one  under  environmental  conditions  corresponding  to  this  particular  oxygen  chemical 
potential.  Using  eq.  (18)  in  section  2.2.1  allows  to  relate  specific  (Tp)-conditions  to  this  chemical 
potential,  and  this  information  can  e.g.  be  included  in  graphs  like  Fig.  4a  in  form  of  additional  x-axes, 
which  give  the  pressure  dependence  at  some  specific  temperature.  Alternatively,  one  can  concentrate  only 
on  the  most  stable  structures,  convert  the  range  of  chemical  potential  in  which  each  is  most  stable  into 
corresponding  (T,p)  ranges,  and  plot  these  stability  ranges  in  surface  phase  diagrams  of  the  form  of  Fig. 
4b.  It  is  important  to  realize  that  both  kinds  of  plots  are  based  on  exactly  the  same  information.  Surface 
phase  diagrams  (Fig.  4b)  provide  a  more  direct  insight  to  the  experimentally  accessible  (T,p)  conditions, 
whereas  free  energy  plots  (Fig.  4a)  summarize  the  two-dimensional  dependence  conveniently  in  the  one¬ 
dimensional,  but  less  intuitive  dependence  on  the  chemical  potential.  Additionally,  it  is  only  in  the  latter 
kind  of  plots  that  also  information  about  the  energetic  difference  to  alternative,  less  stable  surface 
structural  models  is  provided.  These  plots  make  also  immediately  apparent  that  the  transition  from  one 
stable  phase  to  another  occurs  within  the  present  framework  always  at  a  specific  value  of  the  oxygen 
chemical  potential,  cf.  Fig.  4a,  which  is  the  reason  why  the  phase  boundaries  in  (T,p)  surface  phase 
diagrams  of  the  type  of  Fig.  4b  exhibit  similar  curvatures  (lines  of  constant  A ju0)  ■  We  note  in  passing  that 
a  third  and  equally  equivalent  way  of  plotting  the  results  would  be  to  plot  the  stability  ranges  in  (7/7)  p)- 
figures,  in  which  case  these  boundaries  between  the  stable  phases  would  then  result  as  straight  lines  [31]. 

Realizing  that  the  definition  for  the  average  binding  energy  at  T  =  0  K  is 


(27)  Eb 


N0\ 


7t0tal (N0 ,Nm)~  7t0tal (0, Nm )  -  ^7‘°2tal 


with  Eb  >  0  for  exothermicity,  we  arrive  at  the  expression 
(28)  AGad(T,p) «  ^Eh  +  (AM  -  N'jE^  +  ^A/r0(2» 


which  has  a  rather  intuitive  structure:  Forming  the  oxidized  surface  by  accommodating  No  oxygen  atoms 
yields  an  energy  gain  of  NoEb  (per  surface  area),  that  is  opposed  by  the  cost  of  taking  these  O  atoms  out  of 
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the  gas  phase  reservoir,  hence  N0AjUo.  The  equivalent  term  (/V  M  -  /V  V1  )E  ‘‘’tal  comes  only  into  play  for 
oxidized  surfaces  where  the  total  number  of  metal  atoms  is  different  to  the  one  of  the  reference  clean 
metal  surface,  and  represents  then  the  cost  of  transferring  the  corresponding  number  of  metal  atoms  to  or 
from  the  reservoir  represented  by  the  metal  bulk. 

On  the  basis  of  eq.  (28)  the  general  structure  of  a  surface  free  energy  plot  for  a  metal  surface  in  contact 
with  an  oxygen  environment  can  nicely  be  discussed.  In  the  limit  of  an  infinitely  dilute  gas  (A jUo  — ►  -go), 
any  surface  structure  containing  oxygen  (A0  ^  0)  will  exhibit  an  infinitely  negative  A Gad,  reflecting  that  it 
is  very  unfavorable  to  maintain  oxygen  adsorbed  at  the  surface  under  such  conditions.  As  intuitively  clear, 
the  clean  surface  will  therefore  always  result  as  most  stable  in  such  environments.  With  increasing  oxygen 
content  in  the  gas  phase,  A y0  will  become  less  negative  and  so  will  the  A Gad  of  oxygen  containing  surface 
structural  models.  Eventually,  one  of  them  will  exhibit  a  A Gad  >  0  and  will  then  become  more  stable  than 
the  clean  surface.  At  corresponding  oxygen  pressures  and  temperatures,  oxygen  is  getting  stabilized  at  the 
surface,  cf.  Fig.  4a,  and  the  governing  factors  for  this  are  immediately  revealed  by  the  structure  of  eq.  (28): 
The  slope  of  each  AGad-line  is  determined  by  No/A,  i.e.  the  more  oxygen  is  contained  in  the  structure,  the 
faster  this  structure  becomes  more  favorable  with  increasing  oxygen  chemical  potential.  For  surfaces 
preserving  the  number  of  metal  atoms  (Am  =A’m),  the  x-axis  intercept  (i.e.  the  moment  when  the  structure 
becomes  more  stable  than  the  clean  surface)  is  reached  at  A ju0  =  -Eb.  A  more  stable  binding  of  oxygen  in 
the  surface  structural  model  will  correspondingly  shift  the  AGad-line  in  the  free  energy  plot  down  and  the 
x-axis  intercept  to  the  left,  and  will  thereby  render  the  structure  more  stable  than  the  clean  surface  at 
already  lower  oxygen  contents  in  the  gas  phase.  A  large  AGad  at  increasing  oxygen  chemical  potential  and 
therewith  the  chance  to  become  the  most  stable  structure  can  therefore  be  reached  by  surface  structural 
models  that  either  offer  a  strong  binding  of  their  oxygen  species  or  contain  a  large  number  of  oxygen 
atoms  per  surface  area.  This  way,  a  surface  structure  that  strongly  binds  a  few  oxygen  atoms  could  for 
example  become  more  favorable  than  the  clean  surface  at  low  chemical  potentials,  while  another  surface 
structure  with  weaker  binding,  but  higher  oxygen  coverage  will  eventually  become  more  stable  at 
somewhat  higher  chemical  potentials  due  to  its  steeper  slope,  cf.  Fig.  4a. 

The  highest  number  of  oxygen  atoms  per  surface  area  is  ultimately  reached  by  bulk  oxide  structures,  i.e. 
when  the  oxygen  content  in  the  environment  is  high  enough  to  create  an  infinitely  thick  bulk  oxide  on  top 
of  the  metal  substrate.  Since  then  A0  — >  go,  the  corresponding  AGad-line  in  the  free  energy  plot  is  vertical, 
cf.  eq.  (28).  The  intercept  of  this  line  with  the  x-axis  is  given  by  the  condition  that  the  bulk  oxide  becomes 
thermodynamically  more  favorable  than  an  equivalent  amount  of  bulk  metal  and  gas  phase  oxygen 


(29)  gM  O  <xgM+yyQ 


where  gM  0  is  the  Gibbs  free  energy  per  formula  unit  of  the  oxide  bulk.  Using  eq.  (17),  this  yields 


(30)  A/r0>i 

y 


&  Mxoy  (r,g)-xgM(r,g)-^‘°2tall  =  -A Hf{T  =  OK) 
y  2  )  y 


where  AHf(T=0  K)  is  the  heat  of  formation  of  the  bulk  oxide  at  7=0  K  [13,17].  For  any  A ju0  higher  than 
this  limit,  the  bulk  oxide  will  always  be  the  stable  phase,  cf.  Fig.  4a. 

After  this  more  general  discussion  we  proceed  with  the  specific  case  of  a  Pd(100)  surface  in  contact  with 
an  oxygen  atmosphere,  to  illustrate  how  the  direct  screening  approach  to  ab  initio  atomistic 
thermodynamics  works  and  what  it  can  contribute  in  practice.  Typical  for  late  transition  metal  surfaces, 
the  interest  in  this  system  comes  from  the  widespread  technological  use  of  Pd,  for  example  in  the  area  of 
oxidation  catalysis  [32].  Although  this  material  is  known  for  its  propensity  to  form  oxidic  structures  in 
technologically -relevant  high  oxygen  pressure  environments,  the  possible  formation  of  sub-nanometer  thin 
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oxidic  films  (so-called  surface  oxides)  has  only  recently  been  addressed  [17,18,33,35,36].  While 
traditionally  such  films  were  conceived  as  closely  related  thin  versions  of  the  corresponding  (known)  bulk 
oxides,  recent  atomic-scale  characterizations  of  initial  few-atomic-layer  thick  oxide  overlayers  especially 
on  Pd  and  Ag  surfaces  have  revealed  structures  that  have  little  resemblance  to  their  bulk  counterparts, 
and/or  are  influenced  to  a  large  degree  by  a  strong  coupling  to  the  underlying  metal  substrate 
[17,19,20,33-37].  Due  to  this  coupling  and  structures  particularly  suited  for  layered  configurations,  one 
may  expect  the  stability  range  for  such  surface  oxides  to  exceed  that  of  the  hitherto  discussed  bulk  oxides 
[18]. 
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Figure  5:  Computed  free  energy  plot  (a)  and  surface  phase  diagram  (b)  for  O/Pd(100),  following 
the  general  style  of  Fig.  4.  Considered  are  two  ordered  adlayers  with  O  in  the  fourfold  hollow 
sites  (p( 2x2),  %  monolayer  (ML)  coverage,  and  c(2x2),  'A  ML  coverage)  and  the  (V 5  x  V5)R27°-0 
surface  oxide  (0.8  ML  coverage).  Note  the  extended  stability  range  of  the  surface  oxide 
compared  to  the  known  PdO  bulk  oxide.  The  total  energies  (DFT-GGA,  PBE)  used  to  construct 
this  graph  via  eq.  (26)  are  taken  from  refs.  [36,40],  the  surface  unit-cell  area  of  Pd(100)  is  7.8  A2. 


In  the  spirit  of  the  direct  screening  approach  we  therefore  consider  here  the  three  experimentally 
characterized  oxygen-containing  surface  structures  to  date,  namely  two  ordered  adlayers  with  O  in  the 
fourfold  hollow  sites  (p( 2x2),  lA  monolayer  (ML)  coverage,  and  c(2x2),  ‘A  ML  coverage)  and  the  (a/5  X 
V5)i?27°-0  surface  oxide  (0.8  ML  coverage)  [38,39].  The  latter  structure  corresponds  to  a  rumpled,  but 
commensurate  PdO(lOl)  film  with  a  strong  coupling  to  the  underlying  substrate  [35].  Evaluating  the 
calculated  DFT  binding  energies  for  these  three  surface  structures  leads  to  the  results  displayed  in  Fig.  5a 
[35,36,40].  They  nicely  follow  the  more  general  structure  discussed  above:  While  the  clean  surface  is  the 
most  stable  structure  at  the  lowest  oxygen  chemical  potentials,  the  p( 2x2)  structure  exhibits  a  higher  AGad 
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for  Aju0  >  -1.35  eV.  Due  to  its  steeper  slope  (higher  coverage),  the  (V5  x  V5)i?27°-0  surface  oxide 
becomes  even  more  favorable  for  A ju0  >  -1.15  eV,  while  ultimately  in  the  most  oxygen-rich  environments 
(A ju0  >  —0.87  eV)  the  PdO  bulk  oxide  results  as  most  stable  phase.  Converting  this  information  about  the 
stability  ranges  of  these  different  phases  by  means  of  eq.  (18)  into  (7]p)-conditions  leads  to  the  plot  also 
shown  in  Fig.  5b. 

Referring  to  the  more  detailed  original  literature  [17,18,35,36],  we  restrict  our  discussion  of  these  results 
here  to  two  noteworthy  points:  First,  as  motivated  above,  there  is  indeed  a  surprisingly  large  range  of 
(7’,/;)-conditions,  where  the  surface  oxide  structure  represents  the  thermodynamically  most  stable 
structure.  In  corresponding  oxygen  environments,  this  sub-nanometer  thin  film  will  thus  eventually  form 
on  time  scales  set  by  possible  kinetic  limitations,  but  never  grow  thicker.  Due  to  this  finite  thickness,  the 
coupling  at  the  oxide-metal  interface  and  an  atomic  structure  that  can  be  quite  different  to  the  one  of  the 
known  bulk  oxides,  one  might  suspect  new  properties  that  are  distinct  to  those  of  surfaces  of  both  bulk 
metals  and  bulk  oxides,  and  could  thus  be  of  potential  interest  for  applications.  Second,  the  c(2x2) 
structure  is  never  a  most  stable  phase.  This  implies  that  the  frequent  observation  of  this  structure  in  ultra- 
high  vacuum  (UFIV)  experiments  [38,39]  is  a  mere  outcome  of  the  limited  O  supply  offered,  as  well  as  of 
kinetic  barriers  to  the  formation  of  the  surface  oxide,  e.g.  due  to  limitations  in  the  O  penetration  at  the  low 
temperatures  employed  (UHV  experiments  are  typically  performed  by  depositing  a  finite  number  of 
adatoms,  rather  than  by  maintaining  a  given  gas  pressure  [41]). 

Instead  of  further  dwelling  on  the  physics  of  this  system,  let’s  return  now  to  the  methodological 
discussion.  As  already  stated  several  times,  the  validity  of  these  results  is  restricted  by  the  limited  number 
of  surface  structural  models  considered  in  the  direct  screening  approach.  Apart  from  that,  uncertainties  are 
introduced  due  to  the  neglected  vibrational  free  energy  contribution  to  A Gad,  as  well  as  due  to  inaccuracies 
in  the  total  energy  difference  entering  eq.  (26).  Inspecting  the  y-axis  scale  of  the  free  energy  plot  in  Fig. 
5a,  we  see  for  example  that  even  small  changes  in  A Gad  of  the  order  of  ~  ±5  meV/  A2  may  still 
considerably  shift  the  (T,p)  stability  ranges  for  the  various  phases  (due  to  the  shifted  crossing  points  of  the 
various  lines  with  different,  but  similar  slopes).  On  the  other  hand,  the  sequence  of  stable  phases  in 
increasingly  oxygen-rich  environments  (clean  surface,  p{ 2x2)  adlayer,  (f5  X  V5)R27°-0  surface  oxide, 
PdO  bulk  oxide)  is  not  affected  by  such  changes  [40].  Both  the  numerical  uncertainty  in  the  total  energy 
difference  due  to  the  finite  basis  set  employed  in  the  DFT  calculations  and  the  neglected  AFvlb'ad  term  are 
of  this  order  of  magnitude  and  in  the  corresponding  light  the  reported  results  have  to  be  seen.  An  even 
larger  uncertainty  in  the  total  energy  difference  may  result  from  the  approximate  exchange-correlation 
(XC)  functional  in  the  DFT  calculations.  While  in  the  present  example  using  local-density  or  several 
generalized  gradient  functionals  still  led  to  the  same  sequence  of  stable  phases  [40],  extreme  caution  is 
advisable  in  general.  As  always,  the  accuracy  level  is  dictated  by  the  questions  one  wants  to  get  answered. 
If  required,  systematic  improvement  on  the  numerical  uncertainty  in  the  total  energy  difference  and  on  the 
vibrational  free  energy  contribution  is  in  principle  always  possible  (albeit  in  practice  at  possibly  high  or 
prohibitive  computational  cost).  Concerning  the  uncertainty  due  to  the  approximate  XC  functional,  at  least 
one  may  compare  the  results  obtained  with  differently  constructed  functionals.  If  doubts  remain,  a  regional 
XC  correction  or  higher-level  electronic  structure  calculations  may  be  necessary. 

As  a  final  point,  we  briefly  comment  on  the  effect  of  the  neglected  configurational  entropy  contribution. 
As  discussed  in  section  2.2.2  at  sufficiently  low  temperatures  this  term  is  quite  small  and  can  therefore 
only  have  an  effect  when  two  competing  AGad- lines  come  very  close  to  each  other  [15].  This  is  the  case  at 
the  transitions  between  stable  phases,  and  in  fact,  the  deliberately  neglected  configurational  entropy  term 
is  the  reason  why  these  boundaries  are  drawn  abrupt  in  the  surface  phase  diagram  in  Fig.  5  -  even  at  the 
highest  temperatures  shown.  In  reality,  finite  phase  coexistence  regions  should  occur  at  finite 
temperatures,  i.e.  regions  in  which  with  changing  pressure  one  phase  gradually  becomes  populated  and  the 
other  one  depopulated.  With  increasing  temperature,  the  width  of  these  coexistence  regions  around  the 
phase  transitions  increases,  until  eventually  there  are  no  pressures  left  in  which  one  still  finds  the  well- 
ordered  surface  structures  now  displayed  in  Fig.  5.  Only  a  proper  evaluation  of  the  configurational  entropy 


2  - 14 


RTO-EN-AVT-142 


Ab  Initio  Atomistic  Thermodynamics  for  Surfaces:  A  Primer 


term,  e.g.  through  Monte  Carlo  simulations,  can  provide  detailed  insight  into  these  order-disorder 
transitions  and/or  the  phase  coexistence  regions  themselves,  and  corresponding  care  has  to  be  taken  in  the 
interpretation  of  results  at  elevated  temperatures,  when  the  configurational  entropy  term  is  neglected  as  in 
the  direct  screening  approach  to  ab  initio  atomistic  thermodynamics  discussed  here  [6], 


3.0  SUMMARY 

A  predictive  modeling  of  materials  properties  requires  a  consistent  treatment  in  the  wide  hierarchy  of 
scales  from  the  electronic  level  to  macroscopic  lengths  and  times.  The  central  idea  of  ab  initio  atomistic 
thermodynamics  is  to  employ  the  information  on  the  potential  energy  surface  provided  by  modem 
electronic  structure  theories,  in  order  to  calculate  appropriate  thermodynamic  potential  functions.  With  the 
latter,  macroscopic  system  properties  at  finite  temperatures  can  immediately  be  discussed.  At  surfaces, 
such  a  thermodynamic  description  can  be  particularly  useful,  since  it  provides  the  possibility  to  suitably 
divide  the  total  system  into  smaller  subsystems  that  are  mutually  (or  partly)  in  equilibrium  with  each 
other.  This  way,  infinite,  but  homogeneous  parts  of  the  system  like  bulk  or  surrounding  gas  phase  can  be 
efficiently  represented  by  corresponding  reservoirs,  which  e.g.  allows  to  address  surfaces  in  contact  with 
realistic  environments. 

In  this  tutorial  text  we  have  focused  on  a  very  simple  realization  of  this  general  scheme,  namely  the  direct 
screening  approach,  to  determine  the  equilibrium  geometry  and  composition  of  a  solid  surface  in  contact 
with  a  given  environment  at  finite  temperature  and  pressure.  For  the  sake  of  clarity  we  considered  the  case 
of  a  monoatomic  metal  and  an  oxygen  atmosphere,  but  the  conceptual  framework  is  readily  generalized  to 
more  complex  systems,  involving  compounds  like  oxides  or  alloys,  or  environments  containing  multiple 
gas  phase  species.  In  the  direct  screening  approach  one  focuses  on  a  number  of  known  (or  possibly 
relevant)  ordered  surface  structures,  and  directly  compares  which  of  them  turns  out  to  be  most  stable  under 
which  (7’,/;)-conditions,  i.e.  which  of  them  exhibits  the  lowest  surface  free  energy  or  Gibbs  free  energy  of 
adsorption.  This  provides  first  valuable  insight  into  the  structure  and  composition  of  the  surface  in  realistic 
or  technologically  relevant  enviromnents  at  virtually  no  extra  computational  cost  compared  to  the 
underlying  electronic  structure  theory  calculations. 

The  major  limitation  of  the  direct  screening  approach  is  that  its  reliability  is  restricted  to  the  number  of 
considered  configurations,  i.e.  the  existence  of  unanticipated  surface  geometries  or  stoichiometries  cannot 
be  predicted.  This  can  only  be  overcome  by  a  proper  sampling  of  configurational  space,  as  e.g.  provided 
by  modem  statistical  mechanics  methods  like  Monte  Carlo  simulations,  which  then  leads  also  to  a  more 
general  and  systematic  way  of  treating  phase  coexistence  and  order-disorder  transitions.  Last,  not  least, 
one  should  always  keep  in  mind  that  (regardless  of  whether  direct  screening  or  statistical  sampling)  ab 
initio  atomistic  thermodynamics  is  -  as  reflected  by  the  name  -  a  thermodynamic  theory  and  as  such 
describes  systems  that  had  infinitely  long  time  to  fully  equilibrate.  It  provides  no  information  on  what  time 
scale  (with  which  kinetic  hindrance)  this  equilibration  took  place.  For  this,  one  necessarily  needs  to  go 
beyond  a  thermodynamic  description  and  explicitly  follow  the  kinetics  of  the  system  over  time. 
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